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Subgrid  scale  modeling  in  solar  convection 
simulations  using  the  ASH  code 

By  Y.-N.  Young,  M.  Mieschf  and  N.  N.  Mansour 


1.  Motivation  and  objectives 

The  turbulent  solar  convection  zone  has  remained  one  of  the  most  challenging  and 
important  subjects  in  physics.  Understanding  the  complex  dynamics  in  the  solar  con- 
vection zone  is  crucial  for  gaining  insight  into  the  solar  dynamo  problem.  Many  solar 
observatories  have  generated  revealing  data  with  great  details  of  large  scale  motions  in 
the  solar  convection  zone.  For  example,  a strong  differential  rotation  is  observed:  the 
angular  rotation  is  observed  to  be  faster  at  the  equator  than  near  the  poles  not  only 
near  the  solar  surface,  but  also  deep  in  the  convection  zone.  On  the  other  hand,  due 
to  the  wide  range  of  dynamical  scales  of  turbulence  in  the  solar  convection  zone,  both 
theory  and  simulation  have  limited  success.  Thus,  cutting  edge  solar  models  and  numer- 
ical simulations  of  the  solar  convection  zone  have  focused  more  narrowly  on  a few  key 
features  of  the  solar  convection  zone,  such  as  the  time-averaged  differential  rotation.  For 
example,  Brun  & Toomre  (2002)  report  computational  finding  of  differential  rotation  in 
an  anelastic  model  for  solar  convection.  A critical  shortcoming  in  this  model  is  that  the 
viscous  dissipation  is  based  on  application  of  mixing  length  theory  to  stellar  dynamics 
with  some  ad  hoc  parameter  tuning. 

The  goal  of  our  work  is  to  implement  the  subgrid  scale  model  developed  at  CTR  into 
the  solar  simulation  code  and  examine  how  the  differential  rotation  will  be  affected  as  a 
result.  Specifically,  we  implement  a Smagorinsky-Lilly  subgrid  scale  model  into  the  ASH 
(anelastic  spherical  harmonic)  code  developed  over  the  years  by  various  authors.  Readers 
are  referred  to  (Clune  et  al.  1999)  and  (Miesch  1998)  for  a detailed  description  of  the 
ASH  code.  This  paper  is  organized  as  follows.  In  §2  we  briefly  formulate  the  anelastic 
system  that  describes  the  solar  convection.  In  §3  we  formulate  the  Smagorinsky-Lilly 
subgrid  scale  model  for  unstably  stratified  convection.  We  then  present  some  preliminary 
results  in  §4,  where  we  also  provide  some  conclusions  and  future  directions. 


2.  The  anelastic  equations  of  motions 

In  the  solar  convection  zone,  the  depth  extends  over  many  density  scale  heights  and 
thus  the  effects  of  stratification  need  to  be  captured  despite  the  fact  that  acoustic 
timescales  are  much  shorter  than  the  large-scale  convection  timescales.  An  appropri- 
ate approach  is  to  filter  out  sound  waves  while  maintaining  the  density  stratification, 
namely  the  anelastic  approximation  first  proposed  by  Gough  (1969)  and  later  adapted  to 
the  solar  convection  zone  by  Gilman  & Glatzmaier  (1981).  In  stellar  models,  a fluid  layer 
is  unstable  if  the  layer  is  superadiabatic.  The  degree  of  superadibaticity  is  quantified  as 
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(Gilman  & Glatzmaier  1981) 


where  d is  the  depth  of  the  convection  zone,  Cp  the  specific  heat  capacity  at  constant 
pressure,  s the  entropy,  and  r the  radial  coordinate.  The  entropy  gradient  is  calculated 
at  some  fiducial  level  in  the  convection  zone.  In  the  sun,  the  departure  from  adiabaticity 
in  the  convection  zone  is  extremely  small  (Christensen-Dalsgaard  et  al.  1993), 

£ < 1(T4.  (2.2) 

If  we  assume  that  the  velocity  in  the  horizontal  (transverse  to  gravity)  directions  is  mostly 
driven  by  pressure  gradients,  the  coupling  between  vertical  and  horizontal  velocities  leads 
to  the  following  relationship  between  the  convection  velocity  v and  sound  speed  cs 

M = - ~ « 1,  (2.3) 

Cs 

where  M.  is  the  Mach  number  of  the  flow. 

The  basic  idea  in  the  anelastic  approximation  is  to  expand  the  variables  in  terms  of 
the  small  parameter  e,  and  collect  terms  of  the  zeroth  and  first  orders.  As  is  usually 
found  in  perturbation  theory,  the  zeroth  order  equations  describe  the  basic  state,  which 
remain  stationary  over  timescales  at  which  the  first  order  variables  vary.  The  fundamental 
equations  are  the  compressible,  stratified  Navier-Stokes  equations  in  a rotating  reference 
frame 


P + (v  - V)v^  = -VP  + pg  + 2v  x fi0  + V ■ T>,  (2.4) 

p6%  + P&V  ' Vs  = ~V  ' q + 2vp  (ebeb  “ ^(V  • v)2)  “ P£>  (2-5) 

g + V-(pv)=  0,  (2.6) 

where  s is  the  entropy,  © is  the  temperature,  and  £ is  the  nuclear  energy  generation  rate. 
The  viscous  stress  V is  defined  as 

V = 2pv  (eij  - i(V  • v)SijSj  , (2.7) 

and  q is  the  composite  non-convective  heat  flux,  defined  as 

q = —KpQVs  — KrpCpVQ.  (2.8) 

and  the  various  transport  coefficient  such  as  the  viscosity  v,  thermal  diffusion  k,  and 
radiative  thermal  diffusion  Kr  are  assumed  to  be  functions  of  the  radial  coordinate  only. 

Upon  substituting  the  variables  expanded  in  terms  of  the  small  parameter  e,  at  ze- 
roth order  we  obtain  the  hydrostatic  equation.  Denoting  the  zeroth  order  terms  with  a 

subscript  0,  the  equations  at  the  leading  order  are 


-po9, 


P0  = Tlpo&o, 


(2.9) 

(2.10) 
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where  A is  the  radial  gradient  of  a mean  turbulent  pressure 

A = ([/>0(v  ■ V)v  - V ■ V - 2p0{v  x fi0)]r)t  • (2.1i) 

The  angular  brackets  in  Equation  (2.11)  denote  horizontal  averages,  and  the  square 
brackets  denote  the  radial  component  of  the  enclosed  vector.  We  evaluate  A only  as  we 
update  the  reference  state  (denoted  with  a subscript  0). 

At  the  next  order,  we  obtain  the  anelastic  equations  for  compressible,  stratified  fluids 


dpgv 

dt 


+ po(v  ■ V)v  = -VP  + V ■ V + pg  + 2p0(v  x Cl0)  + A r, 


ds 


/9o©owr  + Po©ov  • V(s0  + s)  = V • (kPoQ0V(so  + s)  + Krp0CpV(& 0 + ©)) 
at 


+2vp0 

V ■ (pov)  = 0, 


(^eijeij  g 


(V-v)1 2  \-poS, 


(2.12) 


(2.13) 

(2.14) 


and 

L = = 1 (2  15) 

po  P0  ©o  'jPo  Cp  K ' 

We  note  that  both  the  reference  state  and  perturbation  terms  have  been  retained  in 
Equation  (2.13)  because  the  gradient  of  the  reference  state  entropy  (so)  varies  in  ampli- 
tude from  0(e)  in  the  convection  zone  to  0(1)  in  the  convectively  stable  region.  This 
scaling  is  also  true  for  the  term  involving  thermal  diffusion  of  the  reference  temperature 
(Krd&0 /dr).  We  retain  the  term  KrdQ/dr  because  we  wish  to  allow  for  small  deviations 
of  the  radiative  heat  flux  from  spherical  symmetry. 

The  viscosity  coefficient  v is  assumed  to  be  a function  of  the  reference  density  po.  The 
functional  dependence  is  obtained  from  mixing  length  theory,  and  the  specific  value  is 
chosen  in  an  ad  hoc  fashion  for  a given  numerical  resolution  and  combination  of  parame- 
ters. Our  main  goal  in  this  work  is  to  replace  these  “mixing  length”  viscosity  and  thermal 
diffusivity  with  those  based  on  Smagorinsky  subgrid  scale  model.  More  details  will  be 
presented  in  §3. 

The  reference  state  is  crucial  in  the  simulation  of  the  anelastic  system.  In  the  ASH  code, 
the  reference  state  is  updated  only  every  few  hundred  steps  of  advancing  the  anelastic 
equations.  During  the  update  process,  the  following  equations  are  solved  simultaneously 
for  the  reference  pressure  Po  and  density  p0: 


1 dso  _ ldlnPo  din  po 
Cp  dr  7 dr  dr  ’ 
dP0 
-dT  = 

Using  the  reference  entropy  gradient  (ds0/dr)  and  gravity  {g(r))  profiles  from  the 
previous  time  step,  Equations  (2.16)  and  (2.17)  are  solved  for  the  reference  pressure 
and  density.  Figure  1 is  an  example  from  one  of  the  ASH  simulations.  With  similar 
reference  states  and  viscosity  profiles  to  those  illustrated  in  figure  1,  simulations  of  the 
ASH  code  have  successfully  generated  differential  rotation  profiles  that  are  similar  to 
the  solar  observations.  An  example  from  ASH  simulations  is  illustrated  in  figure  2.  The 
numerical  resolution  for  figure  2 is  65  x 128  x 256  in  the  radial,  latitudinal,  and  azimuthal 


(2.16) 

(2.17) 
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FIGURE  1.  Reference  state  from  ASH  simulations,  the  horizontal  axis  is  radial  coordinate  in  units 
of  solar  radius.  Panel  (a):  Reference  pressure  (dyne/cm2)  within  the  convection  zone.  Panel  (b): 
Reference  density  (g/cm3).  Panel  (c):  the  “mixing  length”  viscosity  (cm2/sec)  as  a function  of 
r. 


directions.  The  left  panel  is  a filled  contour  plot  of  the  averaged  angular  velocity  as  a 
function  of  latitude  and  radius,  while  the  right  panel  shows  the  angular  velocity  at  five 
latitudes  as  a function  of  radius  from  the  bottom  to  the  top  of  the  convection  zone.  For 
this  particular  case,  the  viscosity  is  a prescribed  function  of  radius  as  shown  in  figure 
1(c),  and  the  Prandtl  number  is  fixed  at  1/4.  This  corresponds  to  the  AB  case  in  Brun 
& Toomre  (2002),  in  which  the  differential  rotation  inside  the  convection  zone  is  most 
prominent  among  all  the  simulation  data  presented  in  Brun  & Toomre  (2002).  Figure 
2 is  from  simulations  conducted  on  the  NASA/Ames  SGI  cluster  with  different  initial 
conditions.  A typical  differential  rotation  from  solar  helioseismic  observations  is  shown  in 
figure  2(c)  (same  as  figure  1(b)  in  Brun  & Toomre  (2002)).  The  close  similarity  between 
observational  solar  differential  rotation  and  the  simulation  data  shows  great  promise  that 
better  results  may  be  obtainable  if  the  conveciton  model  is  refined  by  building  in  more 
physics.  For  example,  the  viscosity  and  thermal  diffusivity  may  be  replaced  by  better 
dynamical  models.  Addition  of  a mechanism  to  couple  the  convection  zone  with  the 
tachocline  below  may  result  in  better  agreement  between  simulation  and  observational 
data.  As  a preliminary  initiative  we  have  undertaken  is  to  first  improve  the  model  for 
viscosity  using  Smagorinsky-Lilly’s  model. 
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Figure  2.  Panels  (a)  and  (b):  Differential  rotation  for  the  AB  case  in  Brun  & Toomre  (2002) 
with  different  initial  conditions.  Panel  (c)  is  the  time-averaged  rotation  rates  from  five  years  of 
GONG  data. 


3.  Smagorinsky-Lilly  subgrid  scale  model 

The  principle  underlying  the  Smagorinsky  subgrid  scale  model  is  the  balance  of  pro- 
duction of  subgrid-scale  turbulent  kinetic  energy  and  dissipation  of  isotropic  turbulence 
energy  at  the  characteristic  eddy  size.  Let  the  subgrid  turbulent  production  rate  be  de- 
noted by  V,  then  for  a given  turbulent  viscosity  vt,  the  subgrid  turbulent  production 
rate  is 

V = 2vTeijeij.  (3.1) 

The  dissipation  at  scale  l is  ~ q3/l,  where  q is  a characteristic  turbulent  velocity.  Finally, 
to  find  an  expression  for  vt  in  terms  of  the  resolved  quantities,  we  utilize  the  Prandtl’s 
assumption:  ur  = Ciql,  and  upon  substituting  q into  the  dissipation  rate  and  equating 
dissipation  and  production  rates,  we  obtain  the  Smagorinsky  subgrid  scale  model 

vt  = (Csl)2y/ 2eijeTj,  (3.2) 


where  Cs  is  a constant  in  the  range  of  0.1  < Cs  < 0.3.  Lilly  (1962)  later  extended  the 
Smagorinsky  model  to  turbulence  in  unstably  stratified  convection.  Due  to  the  stratifi- 
cation an  additional  parameter,  the  Richardson  number,  appears  in  the  expression  for 
the  eddy  viscosity.  The  idea  that  led  to  Lilly’s  eddy  viscosity  is  very  similar  to  the  above 
argument  underlying  the  Smagorinsky  model:  buoyancy  production  or  consumption  must 
also  be  accounted  for  in  the  subgrid  energy  balance.  Thus  in  the  stratified  case,  the  dis- 
sipation consists  of  contribution  from  both  eddy  dissipation  and  thermal  diffusion.  The 
Smagorinsky-Lilly  eddy  viscosity  is  simply  the  modified  Smagorinsky  model  with  the 
inclusion  of  stratification  effect 

vt  = (^sO  \/2eyeij 

Here  k?  is  the  eddy  thermal  diffusivity,  and  the  flux  Richardson  number  is  defined  as 


i Kt  d- 

1 Rif 

l 


(3.3) 


_ 9 d\n6 

Kit  = — . 

2 eijeij  or 


(3.4) 
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Figure  3.  The  left  pannel  is  the  reference  entropy  gradient  (erg/g  K m)  and  the  right  pannel 
is  the  averaged  Smagorinsky-Lilly  viscosity  (cm2 /sec)  as  functions  of  r in  units  of  solar  radius. 


Alternatively,  the  flux  Richardson  number  can  also  be  defined  through  entropy  as 


Ri  f = 


9 

dj 


dso/dr 


(3.5) 


where  dso/dr  is  the  reference  entropy  gradient.  In  spherical  geometry,  we  choose  the 
length  scale  l as 

A(r)r2 

Lmaxi.L-max  + i) 

where  A(r)  is  the  radial  grid  spacing  at  location  r,  Lmax  is  the  maximum  (latitudinal) 
angular  mode  index. 

Figure  3 shows  the  reference  entropy  gradient  and  the  averaged  Smagorinsky-Lilly 
viscosity  as  a function  of  radius.  Near  the  bottom  of  the  convection  zone,  the  reference 
entropy  gradient  is  positive,  implying  a stably  stratified  region.  This  corresponds  to  the 
small  bump  seen  in  the  Smagorinsky-Lilly  viscosity  profile. 

In  the  simulations  we  start  from  a stationary  configuration  with  a negative  reference 
entropy  gradient  throughout  the  convection  zone.  The  boundary  conditions  for  the  ve- 
locity are  stress  free,  which  are  not  necessarily  realistic  but  numerically  convenient.  As 
we  are  mostly  interested  in  the  differential  rotation  profile  in  the  convection  zone,  the 
parameters  and  initial  reference  state  are  chosen  to  be  the  same  as  the  AB  case  in  Brun 
& Toomre  (2002).  We  first  obtain  a turbulent  convection  zone  using  the  mixing  length 
viscosity  and  thermal  diffusivity.  As  the  simulation  progresses,  the  differential  rotation 
profile  reaches  a statistically  stationary  state.  We  then  replace  the  mixing  length  vis- 
cosity and  thermal  diffusivity  with  the  Smagorinsky-Lilly  model  with  some  prescribed 
constant  coefficient  Cs.  At  present  we  average  Smagorinsky-Lilly’s  eddy  viscosity  over 
latitudinal  and  azimuthal  angles.  Strictly  speaking  this  is  inconsistent  with  the  energy 
balance  principle  that  gives  rise  to  the  eddy  viscosity.  We  are  now  implementing  a local, 
three-dimensional  version  of  the  Smagorinsky-Lilly  model,  which  requires  a considerably 
tedious  alteration  of  the  ASH  code.  For  the  following  results,  we  present  ASH  simula- 
tions using  the  averaged  Smagorinsky-Lilly  eddy  viscosity  as  a test  study  to  see  how  the 
differential  rotation  may  differ  from  that  obtained  with  the  “mixing  length”  viscosity 
and  thermal  diffusivity. 
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Figure  4.  Eddy  viscosity  (cm1 2/sec)  using  the  length  scale  defined  in  Equation  (4.1)  (solid  line) 
and  the  mixing  length  viscosity  (dashed  line).  On  the  right  is  the  functional  form  of  f(r)  used 
in  constructing  the  eddy  viscosity. 


4.  Preliminary  results  and  future  directions 

Our  experiences  so  far  show  that  the  averaged  Smagorinsky-Lilly  eddy  viscosity  (as 
shown  in  figure  3)  is  too  small,  and  simulations  (with  resolutions  65  x 128  x 256)  blow  up 
soon  after  the  Smagorinsky-Lilly  viscosity  is  activated.  A significant  difference  between 
the  mixing  length  viscosity  and  the  averaged  Smagorinsky-Lilly  viscosity  is  that  the  eddy 
viscosity  is  much  smaller  at  the  boundaries.  Similar  issues  arise  in  large  eddy  simulations 
of  the  atmospheric  flows,  and  it  is  usually  necessary  to  incorporate  a smooth  transition 
from  the  eddy  viscosity  to  a wall-model  viscosity  near  the  boundary.  In  the  solar  case,  we 
can  combine  the  mixing  length  viscosity  with  the  eddy  viscosity  using  a choice  of  length 
scale  l as  follows 

l2  = f(r)l2SL  + (l-m)l2M,  (4.1) 

where  Isl  is  the  length  scale  in  Equation  (3.6)  and  Im  is  defined  as 


1 ( I'M 

Cs  ^^/2eyejjVl  _ W 


(4.2) 


Here  f(r)  consists  of  two  error  functions  so  that  f(r)  = 1 in  the  bulk  of  the  convection 
zone,  and  /(r)  = 0 at  the  boundaries  of  the  convection  zone.  Figure  4 shows  the  resul- 
tant eddy  viscosity  (left  panel)  for  a given  choice  for  /(r)  (right  panel).  Figure  5 shows 
the  differential  rotation  at  a instant  in  time  from  simulations  using  the  length  scale  in 
Equation  (4.2). 

It  is  clear  from  comparing  figure  5 and  figure  2 that  the  differential  rotation  rates 
are  very  different  near  the  top  of  convection  zone.  In  fact  the  solar  observation  suggests 
that  the  rotation  rates  decrease  near  the  top  of  convection  zone,  which  is  the  opposite 
of  the  trend  in  figure  5.  However,  near  the  bottom  and  in  the  bulk  of  the  convection 
zone,  the  rotation  rates  from  the  ASH  simulations  using  the  combined  eddy  viscosity 
are  very  similar  to  those  in  figure  2.  This  implies  that  we  may  expect  to  find  similar 
differential  rotation  profile  once  we  implement  a better  model  than  Equation  (4.2)  for 
the  eddy  viscosity  near  the  top  of  convection  zone.  Better  results  may  be  obtained  as 
well  once  we  implement  the  fully  local  Smagorinsky-Lilly’s  eddy  viscosity. 
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r/Ro 

Figure  5.  Panels  (a)  and  (b):  Differential  rotation  using  the  eddy  viscosity  with  a transition 
to  the  mixing  length  viscosity  near  the  boundaries  of  convection  zone. 
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